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CnJ I ABSTRACT 

The manifestation of point sources in the upcoming Planck maps is a direct reflection of 

the properties of the pixelized antenna beam shape for each frequency, which is related 

| to the scan strategy, pointing accuracy, noise properties and map-making algorithm. 

In this paper we firstly compare analytically two filters for the Planck point source 

extraction, namely, the adaptive top-hat filter ( ATHF) , the theoretically-optimal filter 

\ (TOF). Our analyses are based on the premise that the experiment parameters the 

■ TOF assume and require are already known: the CMB and noise power spectrum 

and a circular Gaussian beam shape and size, whereas the ATHF does not need any a 

priori knowledge. The analyses show that the TOF is optimal in terms of the gain after 

the parameter inputs. We simulate the Planck HFI 100 GHz channel with elliptical 

C*~) ' beam in rotation to test the efficiency of the TOF and the ATHF. We also apply 

the ATHF on the WMAP Q-band map and the derived map (the foreground-cleaned 

map by Tegmark, de Oliveira Costa & Hamilton) from the WMAP 1-year data. The 

uncertainties on the angular power spectrum will hamper the efficiency of the TOF. 

To tackle the real situations for the Planck point source extraction, most importantly, 

Qh the elliptical beam shape with slow precession and change of the ellipticity ratio due 

JL , to possible mirror degradation effect, the ATHF is computationally efficient and well 

j_i ■ suited for the construction of the Planck Early Release Compact Source Catalogue, 

-i— > ■ 
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»,1 INTRODUCTION which will be released in roughly 7 months when the full sky 

is covered once. 

The observation of cosmic microwave background (CMB) Recently, the NASA Wilkinson Microwave Anisotropy 

radiation by ESA Planck mission will be able to provide Probe {WMAP ) experiment (Bennett et al. 2003a) has 

the understanding of our Universe with many key cosmo- demonstrated the importance of the component separation 

logical parameters with unprecedented precision. The sig- tools for extraction of CMB power spectrum from the raw 

nal received from the sky, however, will contain not only data. One of the sources with localized peculiarities of the 

the CMB radiation, but also foreground emissions such as CMB sky is from the point source contamination. For the 

synchrotron, dust, free-free emission, Sunyaev-Zel'dovich ef- WMAP mission, the catalogue of the radio point sources was 

feet from galaxy clusters and extra-galactic point sources. obtained (Bennett et al. 2003b; Pierpaoli 2003) and the final 

Removing foreground contaminations is therefore one of CMB maps were cleaned from the point source contribution, 

the main tasks for the Planck mission. Of these different The WMAP channels at high frequency correspond to 

foreground emissions, the purpose for extraction of extra- the Planck LFI frequency range, in which the contamination 

galactic point sources has two folds: to separate the fore- of point sources seems 3 times higher than that from the 

grounds from the underlying CMB signals to achieve the sci- WMAP . For the Planck HFI channels, the point sources 

entific goal of the Planck mission as this type of foreground are expected to be more significant regarding the issue of 

contamination is rather non-Gaussian. On the other hand, component separation and of point source catalogue, 

the whole sky coverage and a wide range of the Planck ob- In this paper, we compare the theoretically-optimal fil- 

serving frequencies provide an opportunity for producing the ter (Tegmark & de Oliveira-Costa 1998) (hereafter TOF) 

extra-galactic point source catalogue, most notably, the pro- and the adaptive top-hat filter (Chiang et al. 2002b) (ATHF) 

posed Early Release Compact Source Catalogue (ERCSC), proposed for extra-galactic point source detection. Theoreti- 
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cally speaking, the filter proposed by Tegmark & de Oliveira- 
Costa (1998) is optimal in terms of the gain factor. The 
gain factor is an indicator of the amplitude enhancement of 
point sources relative to the background. By 'theoretically' 
we mean that the required information input for this fil- 
ter is known: the CMB and the noise power spectrum, and 
the circular beam shape and size. This filter is known as a 
'Matched Filter' in the field of signal processing. 

A lot of effort has been devoted to the development of 
other linear filters (Cayon et al. 2000; Vielva et al. 2001; 
Sanz, Herranz & Martmez-Gonzalez 2001; Vio, Tenorio & 
Wamsteker 2002; Barreiro et al. 2003; Vielva et al. 2003), 
namely the spherical Mexican Hat Wavelet (MHW), and 
the so-called pseudo-filter. Vio et al. (2002) and Barreiro et 
al. (2002) have performed the comparison between the TOF 
and the pseudo-filter. Vielva et al. (2001, 2003) have applied 
the MHW on simulation data for the Planck mission. 

The Planck in-flight beam shape properties are cru- 
cial in connection with the point source extraction. There 
are many issues regarding the in-flight antenna beam shape 
properties and its reconstruction (Burigana et al. 2001; Fos- 
alba, Dore & Bouchet 2002; Chiang et al. 2002a; Nasel- 
sky et al. 2002). First of all, according to optics calcu- 
lations the main beams are not circular Gaussian (Burig- 
ana et al. 2001), but elliptical, which is normally formulated 
as bivariate Gaussian functions. Secondly, with the most 
likely scan strategy, the inclination angle of the beam is not 
all parallel through the sky due to possible precession of its 
spin axis, but is rather slowly rotating along circular scans. 
Note that this slow rotation of the antenna beam is differ- 
ent from the balloon-borne experiments such as MAXIMA- 
1, which has a fast rotating beam (Wu et al. 2001), or 
WMAP, in which the pixels have hits of multiple orientations 
(Page et al. 2003) , and therefore a resultant quasi-symmetric 
pixel-beam. Furthermore, there could exist mirror degrada- 
tion effect during the approximately 15-month routine op- 
eration of the Planck mission (Naselsky et al. 2002) , which 
will mimic the change of the inclination angle and the shape 
of the beam. The combined effect from the above-mentioned 
factors will cause the in-flight beam shape to change with 
time not only its ellipticity ratio <t+/<j_, where a+ and cr- 
aze the major and minor axis of the ellipse, but also the 
inclination angle. There are also some minor effects such as 
the issues on pointing accuracy and the noise properties that 
could also mimic the beam degradation. The filters which as- 
sume a circular Gaussian beam for point source extraction 
will need corresponding modifications. 

The adaptive top-hat filter is flexible for real situations. 
The ATHF is similar to the so-called Gabor transform of 
the signal, using a adaptive width of the kernel. The ATHF 
has a top-hat shape with two cut-off scales working in har- 
monic domain. The two cut-off scales ^ max and £min serve to 
cut down the unwanted power from both sides of the power 
spectrum. Simultaneously, the two scales retain the essen- 
tial part of power spectrum where the convolution (by the 
beam) of point sources is most pronounced. 

We also propose a median filtering technique for remov- 
ing the point source from the maps. The median filter can 
be used on combination with any methods of point source 
detection in order to increase the precision of the CMB map 
reconstruction. 

The layout of this paper is as follows. In Section 2 we 



discuss the antenna beam shape properties by firstly intro- 
ducing the Planck scan strategy and the definition of the 
window function. In Section 3 we present the detailed anal- 
ysis on the theoretically-optimal filter (TOF) (the Matched 
Filter) and the adaptive top-hat filter (ATHF). In Section 4 
We compare by simulations the TOF and the ATHF on gain 
capability on more realistic situations and apply the ATHF 
on Discussions and conclusion are in Section 5. 



2 THE PLANCK ANTENNA BEAM SHAPE 
PROPERTIES 

2.1 The scan strategy 

The manifestation of point sources in the upcoming Planck 
maps is a direct reflection of the properties of the pixelized 
antenna beam shape, which is related to the scan strategy, 
pointing accuracy, map-making algorithm and the extrac- 
tion of the systematic effects from the time-ordered data 
(TOD) set and pixelized maps as well. Before our analysis of 
point source problem, in particular, for the Planck Early Re- 
lease Compact Source Catalogue (ERCSC), we need to de- 
scribe the Planck experiment when the TOD contain the in- 
formation about the signal (and noise) from a large number 
of circular time-ordered scans (Delabrouille, Patanchon & 
Audit 2002; Chiang et al. 2002a). Below we assume for sim- 
plicity that the systematic features are already removed and 
the instrumental noise is Gaussian white noise. In the tem- 
poral domain the observed signal mj is the combined sig- 
nal dt of the CMB and foregrounds from the sky (here- 
after beam-convolved 'sky signals'), plus random instrumen- 
tal noise n t , 

m t = d t + n t (1) 
where 

d t = ^ B t ,t m ai m Yi m {r t ), (2) 

1=0 m=-l 

and Bt,tm is the multipole expansion of the time-stream 
beam B t (r t ). In Eq. (2) ae m are the corresponding multipole 
coefficients of the CMB and foreground signal expansion on 
the sphere and Yt m are the spherical harmonics. 

Following Tegmark & Efstathiou (1996) we also assume 
that map making algorithm is linear. Thus the signal in each 
pixel s p should have the following relation 

JV 

dt = ^M iiP s p , (3) 

where Mt, p is the corresponding pointing matrix and s p is 
the sky signals convolved by the pixel beam B p j m , 

oo I 

{**), (4) 

e=o m=-l 

Xp is the two-dimensional vector with its components x p and 
y p denoting the location on the surface of the sphere. The 
definition of the pointing matrix Mt jP depends on the scan 
strategy of an experiment. Below, as a basic model, we will 
use the model of the scan strategy of the Planck mission 
which was discussed by Burigana et al. (1998). Let s be the 
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unit vector along the satellite spin-axis in the anti-Sun di- 
rection and o is that along the direction of the optical axis of 
the telescope. The angle between spin-axis and optical axis 
is o?o = cos _1 (s ■ o) ~ 85°. The satellite itself will scan the 
same circle 60 times around the spin-axis at Q = 1 r.p.m.. 
Each hour the spin-axis is manoeuvred along the ecliptic 
plane by ~ 2.'5. Note that for our geometrical model, during 
the circling of the optical axis around the spin-axis in each 
step, the inclination of the beam relative to the optical axis 
are stable. For modelling of the Planck antenna beam shape 
we will describe the simplest scan strategy with the spin axis 
right on the ecliptic plane without any additional (regular) 
modulation (e.g. it can be above or below the plane due 
to precession). This model reflects the geometrical proper- 
ties of the beam asymmetry and their manifestation in the 
pixelized maps. Under the assumptions mentioned above we 
can now describe the elliptical beam shape model as follows, 



B t (x — x t ) = exp 
where 



-(RU) T D _1 (RU) 



U 



x — x t 

y - yt 



(5) 



(6) 



R is the rotation matrix which describes the inclination of 
the elliptical beam, 



R 



cos a 
— sin a 



sin a 

cos a 



(7) 



with a being the inclination angle between x axis and the 
major axis of the ellipse. The D matrix denotes the beam 
dispersion along the ellipse principal axes, which can be ex- 
pressed as 



D = 









2 



(8) 



The center of the Cartesian coordinate system is denoted by 
Xt and y t at some moment t with z axis along the o, and 
x and y axis on the plane tangent to the celestial sphere. 
We will choose below the standard orientation of the lo- 
cal Cartesian coordinate system with x axis parallel to the 
scanning direction. 

After some time At = 2-KN TOt /fl, where iV ro t = 60 is 
the number of the sub-scans, the spin axis of the satellite 
shall be stepped along the ecliptic plane to s with the small 
deviation As = s — s ~ 2'.5. For a small part of the sky we 
can neglect the rotation of the Cartesian coordinate system 
and choose the x and y axes parallel to x, y axis. 

The received signals in each pixel p depend on the orien- 
tation of the pixel beam B p (x) and the location of its center 
(Wu et al. 2001). For the upcoming Planck non-symmetric 
spatially-dependent beam, the convolution of the signal with 
asymmetric beam immediately produces non-Gaussian sig- 
nature coupled with the underlying signal, which affects the 
estimation of angular power spectrum. Although the real 
beam profile (including the sidelobes) are complicated, we 
begin our analysis from the model of the time-dependent, 
elliptical beam shape which is applicable up to ~ —30 dB 
level. At low amplitude part of the beam shape (< —30 dB) 
the complexity of the beam estimation increase dramati- 
cally. Moreover, due to possible precession of the spin axis s 
(around the standard orientation), which can be described 



as some additional noise, it is possible to find some trans- 
formation of the width of the beam taking into account the 
statistical properties of the precession. In such a case the 
angle between the spin axis and the optical axis should be 
4>(t) — 4>o + v(t), where 4>a ~ 85° = const, and v(t) <C <j>o 
denote a random precession. 

2.2 The window function 

Due to the scan strategy and (possible) precession, the pixel 
beam pattern, in principle, should be different for different 
pixels on the sphere, which means that we need some mod- 
ification of the standard tools for the Ce power spectrum 
extraction by taking into account any possible peculiarities 
from the bivariate Gaussian beam shape. The finite resolu- 
tion of the pixel beam means that the signal we receive from 
the sky is smoothed by the beam, 



T(p) 



I 



dQ n B(p, n)TcMB(n) 



(9) 



Usually the influence of the pixel beam on the observed 
signals can be described in terms of the window function 
(White & Srednicki 1995; Souradeep & Ratra 2001) 



w; 



47T 



(P.<Z) 



2£ + 



Bim,pB} mt q 



(10) 



where p and q denote the positions of the corresponding 
pixels in the map and 



B e , 



,p = J 



dttnB(p,n)Y e * m (n) 



(11) 



is the spherical harmonic transform of the beam function 
pointing in the direction p. This window function determines 
the CMB signal correlation matrix as 



(P,Q) 



oo 

_ ^ 21 + 1 a 



(12) 



In the framework of the Planck mission it is possible to use 
the flat sky approximation (as a part of the whole sky Ce 
harmonic analysis) in order to estimate the point source con- 
tamination and beam influence on the high multipole part 
of the Ce power spectrum. For such flat sky limit the com- 
putation cost is minimal, as we can implement Fast Fourier 
Transforms (FFT) for the CMB anisotropy. When the sky 
signal is convolved with an elliptical beam with an angle be- 
tween its major axis and the scan direction, the temperature 
becomes 

T{P) = J™ J^y^ k ' PB ( R a [k])TcMB(k), (13) 

where B is the Fourier transform of the beam profile at the 
origin, and R a [k] is the rotation operator with an angle a 
with respect to the scan direction (Souradeep & Ratra 2001), 
which can be written explicitly as 



B(R a [k]) = exp 



-i(RK) T D(RK) 



where R and D are defined as in Eq. (7) and (8) and 



(14) 



(15) 



© 0000 RAS, MNRAS 000, 000-000 



4 Lung- Yih Chiang & Pavel Naselsky 



The window function now can be expressed as (Souradeep & 
Ratra 2001) 



Wfc(a,A0) = 

f -2n , , 



2% 



k e ikAe B(R a [k])B*(R a+Ae [k] 



(16) 



where 4>k is the angle of wave vector k = (k x ,k y ) in the k 
space such that k x = kcos<j> k and k y = fcsin0 fc . Moreover, 
for the flat sky approximation (and without the influence 
of the instrumental noise and any peculiarities of the scan- 
ning such as precession) we can use the assumption of sta- 
ble beam orientation, a — const., with the corresponding 
Fourier image for the elliptical part of the beam shape, 



B(k) 



exp 



o\ (k x cos a + k y sin a) 2 



a'- {—k x sin a + k y cos a) 2 



— exp 



kl (ky-k x b ^fa 2 a 4 



2a? 



(17) 



where a 2 — a+a-,r = a+/a-,a 2 = (cos 2 a + r 2 sin 2 a) /ra 2 
and b 2 — — (r 2 — 1) sin2a/2r<r 2 . Without loss of generality, 
we can rotate the coordinate system such that a = 0, and, 
for small AO, AO = = (9 x ,9 y ) with f3 = axctan(9 y /9 x ). 
The window function can be written as 

W 2 k {AO) ee W 2 {6,13) = /; cos[fc0 cos(0 fe - (3)] 

x cxp{-fc 2 cr 2 r [l + (r- 2 -l)sin 2 (? !> fe ]}, (18) 

with the bivariate beam profile 



B(k) = exp [1 + (r- 2 - 1) sin 2 cfk] } 



(19) 



The integral in Eq. (18) can be expressed analytically 
(Souradeep & Ratra 2001). If the asymmetry of the beam 
shape is low then we can express Wk (9) as an infinite series 
in powers of asymmetry parameter tfe = k 2 a 2 (r~ 2 — 1), 



the influence of the anisotropy of the beam on the CMB sig- 
nal after convolution which transforms isotropic CMB signal 
to anisotropic one. 

For applications of the above-mentioned properties of 
the beam shape and the window function for point source 
extraction, let us define the mean beam profile as follows, 



2t + 



(23) 



In the flat sky approximation, taking Eq. (19) into account, 
the mean beam profile is related to the ellipticity parameter 
r as 



B(k) = exp 



fcvy_+i) 

4r 



lo 



k 2 a 2 {r 2 - 1) 



4r 



(24) 



where Io(x) is the Bessel function of the first kind. For the 
window function W 2 {k) in the flat sky approximation we 
have 



W 2 {k) 



r 44 

Jo 11 



— exp 



exp 



k 2 a 2 



[l + (r 2 - 1) sin 2 <(>} 



k 2 o 2 (r 2 + 1) 



2r 



feV(r 2 _ i) 



2r 



(25) 



The last two approximations play a crucial role in point 
source extraction. First of all, the actual value of the beam 
profile at the location of the point source is related to the 
orientation of the beam as 



B(k) = B(k)Io 



k 2 a 2 {r 2 - 1) 



exp 



4r 

k 2 o 2 (r 2 - 1) 
4r 



cos 20fc 



(26) 



and is different from the mean beam profile. These differ- 
ences depend on the position of the point source in the map 
and for a given value of the flux Si the result of the point 
source subtraction should be different for different points of 
location. 



w^e,/3) = J2i-^n^) n w k , 



where 



(20) 



iffc.n = B Q,n + (cos/3) 2n exp(-fe 2 o- 2 r) 
" (tan P) 2m F 2 (± + m, ±, 1 + n; -k 2 9 2 /4) 



m=0 



r(n - m + l)r(m + 1) 



(21) 



r(a) and B(a, b) are the Eulcr Gamma function of the first 
and the second kind, respectively, and F2 is a generalized 
hyper-geometric function. For beams with symmetric shape 
where r = 1 and er+ = u_ = a, we obtain from Eq. (21) 



Wl(6) = J (k9)e 



(22) 



where the zeroth-order Bessel function Jo(k6) corresponds 
to the standard asymptotic of the Legendre polynomials 
P e (cos9) w Jo(k9) at k = £ + 1/2. In a general case when 
r > 1 the window function W 2 (9) is a function of the separa- 
tion angle 9 and phase (3 as from Eq. (21). This fact reflects 



3 OPTIMIZATION OF THE LINEAR FILTERS 

In this section we would like to elaborate for point 
source extraction the subtleties of the 3 linear filters: the 
theoretically-optimal filter, the adaptive top-hat filter and 
the pseudo-filter. The analyses and comparisons are based 
on the assumption that the beam shapes for all frequency 
channels are circular Gaussian and are known if required 
(by the TOF and the pseudo-filter) and the signal proper- 
ties, both the CMB and pixel noise power spectrum, arc 
known if required (by the TOF). The analysis and compari- 
son is in terms of the gain factor. The gain factor is defined 
as Rt/Ri, where Rf is the signal Sf to noise Nf ratio with 
"f" indicating after convolution by the filter. 



3.1 The theoretically-optimal filter 

In the following analyses for the theoretically-optimal filter 
(TOF) we adopt the convention of notation by Tegmark 
& Efstathiou (1996), Tegmark & de Oliveira-Costa (1998) 
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(hereafter TD98), i.e., the de-convolved power spectrum of 
the combined signal is denned as 



Qsky 



c p 



(27) 



where C| ky is the sky signals including CMB and all fore- 
grounds but the point source contribution, C plx is the noise 
power spectrum, which is taken as (FWHM<j p j x ) 2 , and Wf 
is the window function. Thus the combined (except point 
sources) power spectrum in the map is C^ ot ~ ps = CtWj . 
Note that the in-flight beam Be m and the window function 
Wf definitely have some error bars caused by possible mirror 
degradation effects, the galactic foreground contaminations 
and the instrumental noise. Thus simple inversion of Wf 2 
for the power spectrum requires specific renormalization due 
to the antenna beam shape properties. 

The idea to construct the theoretically-optimal filter 
(TOF) for point source extraction is outlined in TD98. We 
can write down the signal as the sum of the point source 
contribution and the rest (TD98) 



x ( r ) = 9^^Si5 D {r t ,r) + V' ai m Yim(r), 



(28) 



where 8d is a Dirac delta function, Si is the flux of the point 
source at the direction n, at m is the coefficients of the signal 
decomposition with (|a^ m | 2 ) = Ce as in Eq. (27), Ye m (r) the 
spherical harmonics, and 



- _L ( hc V 

9 ~ 2k V fcTcMB / 



[2 sinh(/i^/2fcT C MB) 
(hv/kTcMB) 4 



(29) 



This filter is named as 'theoretically-optimal' because it is 
constructed in such a way that the gain factor will be max- 
imal, i.e., when the 'observed' signal, B ® x(r) is convolved 
with the filter F, 



V{r) = g^Si(F®B)(n-r) 



only ^-dependence, but not on the azimuthal numbers m. 
This condition is only valid when the sky signal and the pixel 
noise are Gaussian. If the sky signal (or the pixel noise) is 
non-Gaussian, however, then Eq. (31) needs the correspond- 
ing modification. For the following analyses we assume for 
simplicity that the Gaussian assumption is appropriate for 
the CMB plus foreground. 

The determination of the TOF is similar to the defini- 
tion of functional derivatives definition. Suppose that F e cor- 
responds to the maximization of the ratio R from Eq. (31), 
then a small variation fe of the shape from the Ft shall shift 
the function R{Fi + fe) away from the maximal value 7? max , 
which is proportional to the fe, fe---, if I f\ *C 1- So we 
have 



F e = F t + ft, 

and obtain the following formulae from Eq. (31), 



A ee R(F e ) - R(F e 



where 



>5> 



M + Q 
aP ' 



M = 2(3 " PW'CnFn] fn, 



21+ 1 
4tt 



B e F e , 



(32) 
(33) 

(34) 
(35) 
(36) 



Q 



E ^r Bnfn 



^ 2 E 



—. — W n f n C n 

4ty 



(37) 



tm 



Fe m Be m aem.Ye m (r), 



(30) 



the maximization of the 'signal' to 'noise' ratio (Sf/Nf) will 
give us the filter shape, where Sf is the square root of the 
variance of the point source amplitudes and Nf now is the 
rms of the combined signal, both convolved by the filter F. 
* Equivalently, we want to maximize 

9 2 (Ei E, SiSjftF ® B)(n ■ r)][(F ® B)( rj ■ r)]> 



R = 



Z e (2l + l)/(^)F?W?Ct 

[Ee W + l)/(^)FeB7\ 2 
J2 e (2£+l)/(4ir)F?W 2 C e ' 



(31) 



where ® denotes convolution, Be m are the coefficients of the 
spherical harmonic decomposition of the beam profile, Fe m is 
the multipole expansion of F , and Be is the mean beam. We 
thus can obtain the TOF for each Planck frequency channel. 

Note that in reaching Eq. (31) the filter is assumed 
isotropic in the Fourier rings and statistically homogeneous. 
Thus in the spherical harmonic decomposition this filter has 

* Note that this definition of our signal-to-noise ratio is inverse 
of that in TD98. 



and 



E 



2£ + 1 
4tt 



W;(F e + fefCe. 



(38) 



Equivalently, to maximize the gain factor R(Fe) is to 
take M — 0, the condition that the first functional derivative 
is zero, 8R/SFe — 0. If we assume that the TOF profile is 
an analytic function, then from Eq. (34) we get 



aB n 



(39) 



As one can see from Eq. (35) and (36) the coefficients a 
and P are also the functionals of the TOF profile. Substitut- 
ing Eq. (39) into Eq. (35) and (36) we obtain the following 
relationship between a and /3 functionals, 



^ 2 = "E 



21 + 1 Be 
4tt WK 



(40) 



Taking into account the definitions of the a and f3 function- 
als and Eq. (40), we reach the conclusion that the TOF F„ 
from Eq. (39) have the form 



pBe 
WjCt ' 



(41) 
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where p — a/@ and (0 1 /a)g 2 ^\ Sf is the maximal gain 
after filtering, 



21 + 1 B e 
4tt WfCt' 



(42) 



If the TOF Fe corresponds to the global maximum of the 
gain (in the class of the analytic functions) then the second 
functional derivative 8 2 R/8F 2 should be negative at F e = 
F(. For small variations around Fe this corresponds to the 
condition Q < 0, according to Eq. (33). After substituting 
fe = Fe ee into Eq. (37) we get 

Q = a ^ ^ h e h n e e e n — ^ h n ^ h e e 2 ) 

\ e n n I J 

= -^S^«(^-e„) 2 , (43) 

"in 

where he = (21+ l)/(47r)£^ 2 /WfCe- As one can see from 
Eq. (43), for any values of the ee functions, the function Q 
is negative and the filter Fe = Fe corresponds to the maximal 
gain in the functional space of the linear filters. 

The filter shape of Eq. (41) for point source extraction is 
a generalization of the TD98 filter, which is obtained under 
the assumption of circular Gaussian antenna beam shapes. 
The new element which we get for the more complicated 
Planck antenna beam shape in Eq. (41) shows that for the 
TOF construction, instead of the TD98 model, we need to 
know the mean beam (over all m from Eq. (23) beam pro- 
file), the window function W 2 , and all the signal properties 
as well. 

One of the most crucial factors affecting the efficiency of 
the TOF is from the systematic effects, namely, the possible 
degradation effect of the primary and secondary mirror sur- 
face during the flight (Naselsky et al. 2002) . As mentioned in 
the Introduction, the mirror degradation could change the 
ellipticity ratio and the beam size, so the calibration and 
reconstruction of the Planck in-flight antenna beam shape 
will have direct consequence of Eq. (23) and the window 
function. 

From a practical point of view, when we apply the TOF 
of Eq. (41) for the Planck point source extraction, in par- 
ticular for the ERCSC, it is necessary to know precisely 
the power spectrum of the sky signals and the pixel noise 
properties. It is, however, only possible at the earliest af- 
ter the analyses of component separation and Ce extraction 
when the mission is completed. Conversely, for the realiza- 
tion of the Ce extraction it is necessary to know the window 
function W e 2 for all frequency ranges. ^ This means that for 
the ERCSC construction it is better to to use primitive fil- 
ters which require as little information as possible about the 
properties of the antenna beam shape and the sky signals. 

One additional possibility for the new class of the lin- 
ear filters comes from Eq. (34) at M = but for the non- 
analytical shape of the optimal filter Fe, which is not in the 
scope of this paper. In the next section we will demonstrate 

t One plausible idea is to use very bright point sources with flux 
above 1 Jy for determination of the window function and the mean 
beam shape at any time of the observations (see, e.g., Chiang 
et al. 2002a). 



that the maximization of the gain factor can be obtained 
using a simple top-hat filter. 



3.2 The Adaptive Top-Hat filter (ATHF) 

Chiang et al. (2002b) introduce another class of filter simi- 
lar to the Gabor transform (Gabor 1946; Hansen, Gorski & 
Hivon 2003): the adaptive top-hat filter (ATHF), which is 
implemented in the harmonic space with two adaptive cut- 
off parameters, F™ — 8(£ — £ m in)6(^max —£), where O is a 
Heaviside function, £ m i n and i? max are two adaptive parame- 
ters. The ATHF allows us to investigate a new class of non- 
analytical functions for the ERCSC. Point source extraction 
by the ATHF is based on the idea that each point source in 
the map (and the TOD as well) manifests itself as a non- 
Gaussian feature which is convolved with the beam shape. 
The goal, therefore, is to distinguish these non-Gaussian fea- 
tures from the CMB signal and pixel noise by using some 
general properties of the CMB map only. The ATHF is gen- 
eralized from the amplitude-phase analysis method for point 
source extraction in the Planck maps (Naselsky, Novikov & 
Silk 2002). In filtering, the parameter i? m i n suppresses all 
the power of the signal from the sky at the multipole range 
i < imin for which the beam shape influence on the signal 
is not crucial. This part of filtering by the ATHF reflects 
the fact that the pure CMB signal has £~ 2 tail, whereas 
dust, synchrotron and free-free emission have ~ l~ z tails for 
the foregrounds (Tegmark & Efstathiou 1996). The param- 
eter £max, on the other hand, suppresses the high multipole 
power from pixel noise and secondary anisotropies, in order 
to relatively increase the point source amplitudes above the 
threshold of detection. 

Non-symmetric beam break the statistical isotropy of 
the Fourier ring in the power spectrum and the ellipticity of 
the beam will manifest itself in the Fourier domain (Chiang 
et al. 2002a). The usual estimation of the ellipticity ratio of 
Planck antenna beam (~ 1.2 -j- 1.3) is well above the ratio 
of the two proposed cut-off scales of the ATHF for all the 
Planck channels (Chiang et al. 2002b) , which means the cut- 
off scales of the ATHF will retain the part of the elliptic 
shape on the Fourier ring. 

Before our analyses, we would like to emphasize again 
the adaptiveness of the ATHF which does not at all need 
any a priori information of the experiment parameters. In 
principle, we can start from any values of the £ m in and £ max 
parameters for investigation of the peak amplitudes in the 
CMB map without any assumptions on the properties of the 
signal and noise. We can choose, for example, £1, = £ p i x , 
where £ p i x is the multipole mode corresponding to the pixel 
size, and £^ in = 2, the lowest multipole mode. After the first 
step of the iteration, the maxima might contain both CMB 
peaks and point sources above some specified threshold vtO~f, 
where 07 is the square root of the variance of the filtered 
map and Ut is a dimensionless amplitude, which, following 
the standard criterion (see, e.g. TD98), we will set vt = 5. 

The next few steps of the iteration are to increase the 
£ m in parameter. Because the power of the low multipole part 
of the CMB and foreground signal are gradually subtracted 
by the increasing £ m in, the new variance of the signal after 
filtering 0% must be smaller than the preceding one, but the 
point source amplitudes changes only slightly, which results 
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in pushing the gain factor higher. The steps of increasing 
£ m in are carried on up to the step where the gain factor dips 
after reaching a maximum value. We can then tune the ! maJ < 
parameter in the similar way by lowering it from l^ix = ^pix- 
Tuning of the cut-off scales of the ATHF will produce 
the so-called Brownian motion of the CMB peaks, as com- 
pared to motions within the beam size for beam-convolved 
point sources (Chiang & Naselsky 2003), which is useful in 
lowering the threshold of the criterion for point source de- 
tection. 

For practical implementation of the iteration scheme, 
however, it is not necessary to start from (£^ in ,£mlx) = 
(2, £pi x ). In order to minimize the time needed for the re- 
alization of the iteration scheme, we can use as the first step 
from the suggested sets of € m in and £ max for all the Planck 
channels (Chiang et al. 2002b). Below to obtain the maximal 
gain from the ATHF, we can give the same treatment as that 
for the TD98 filter, but treating fc m in and fcmax as variables. 
Again, the same assumptions are made on the properties 
of the signal and noise similar to the TD98 filter and its 
generalization TOF [Eq. (39)], namely, the mean beam, the 
window function and the power spectra of the signal and 
noises are known. For analysis we will use the flat sky ap- 
proximation for a small path of the sky. The generalization 
of the analysis for the whole sphere can be done without 
special modification. 

Let us start with Eq. (31) for the flat path of the sky, 
the cut-off shape of the ATHF transforms the fcmin and fcmax 
parameters on to the summations at both numerator and 
denominator, 



HFI 100 GHz 



c, 



c pix /w t 2 




1000 

I f 

min max 



Figure 1. The simulated angular power spectrum of the CMB 
and pixel noise for the Planck HFI 100 GHz channel. The solid 
curve is the power spectrum C| ky , which includes the CMB signal 
and foregrounds, and the dash line is the dc- convolved noise power 
spectrum C p '*/Wf. The thick curve represents the de-convolved 
power spectrum of the combined signal Ce ■ The shaded area shows 
the optimal top-hat filtering range (^ m in,^max) for this channel, 
taken from Chiang et al. (2002b). We can see from the figure that 
C,(W) ~ C7P ix /IT 2 (W) and C*(* min ) ~ C sk ^ min ). 



-Rf^min? ^max) — Q ^ ^ S{ 



f dkkB(k) 



^max 

/ dkk W 2 (k)C(k) 



(44) 



and 



I .till N 

) / dkkB{k), 

^min 



(46) 



where the shape of the filter F(k) — Q(k — fc m i n )0(fc max — fc) 
now is quite peculiar. By the definition, the perturbations 
of the fcmin and fcmax parameters, <5fc m m and 5fc m ax, produce 
the perturbations of the filter shape 



/(fc) — @(fc fcmax T <5fc m ax)@ (fcmax fc) 

+ e(fc-fc m l n + Sk m l n Mfcmin-fc) 



(45) 



As one can see from Eq. (45) these perturbations corre- 
spond to \ \,f(l) | 3> \F(l)\ | at the ranges fc ma x — <5 fcmax < k < 
fcmax and fcmin — <5fc m in < fc < fc m in- Thus finding the optimal 
values of the fcmin and fcmax parameters through the same 
treatment as in the last section is not appropriate. From a 
theoretical point of view that means that for a given top-hat 
shape of the filter the point of maxima of the gain parameter 
does not correspond to any solutions for Eq. (34) . 

In order to find the optimal fc m in and fcmax values we use 
instead the following conditions d_R(fc m in, fcmax)/<9fc m in — 
and 9i?(fc m in, fcmax)/9fc m ax = 0. For Eq. (44) these condi- 
tions lead to 



"rn; 

/ 



2£(fcmax) / dkkW 2 (k)C{k) 



2B(fcmin) J dkk W 2 (k)C(k) 

femin 

^max 

= W 2 (fcmin)C*(fcmin) j dkkB{k). (47) 
&min 

We get from Eq. (46) and Eq. (47) 

IF 2 (fcmin)C(fcmin) = E.^™'™'. \ w 2 (fc ma x ) C(fc m ax ) ■ (48) 
B (fcmax) 

Apart from a trivial solution in Eq. (48), with undefined 
value of the fc = fc max = fcmin, there are non-trivial 
solutions. For comparison with other filters we use the 
same circular Gaussian model for the antenna beam and 
the window function, W 2 (k) = "B (fc) = exp(-fc 2 6»|), 
where B is the FWHM/ v / 81n2, so from Eq. (48) we have 

-B(fcmax)C*(fcmax) = B (fcmin) C^(fcmin) ■ We can assume that 
fcmin^B < 1, 

B(fcmin) = exp(-fc mi X/2) ~ 1 » exp(-fcLx#!/2). (49) 

We further have the following two assumptions. At the 
range of fc <J fc m in, the power spectrum of the combined 
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Figure 2. The point source extraction rate from the TOF and 
the ATHF. We put 50 point sources with 1.5(tcmb randomly 
in the CMB map simulating the Planck HFI 100 GHz channel 
with the corresponding pixel noise. This number of point sources 
is too large for the 100 GHz frequency channel, it nevertheless 
serves as the efficiency (percentage) of extraction rate. The el- 
liptical beam is rotating across the simulation patch with cllip- 
ticity ratio 9+/ 9— = 1.3. The plot is the extraction percentage 
of point sources against the 9 inserted to the TOF (in terms of 
00 = 10.7/ V8 In 2 (arcmin) ). In the dotted curve, the power 
spectrum is precisely known. The short-dash curve is the case 
when the slope of the CMB power spectrum has 5 per cent more 
tilt, whereas the long-dash line 5 per cent less tilt. The horizon- 
tal level marks the ATHF extraction percentage. With the exact 
power spectrum, the TOF is still able to detect 85 per cent of the 
point sources if the supplied 9 is chosen in between 9+ and 9- . 



signal C(fcmin) is mainly determined by the sky signals, and 
the pixel noise contribution is not significant, so 



C(k n 



C sky (k n 



(50) 



On the other hand, for the range of k ^ fc max , the power 
spectrum is mainly determined by the pixel noise power 
spectrum, which leads 



C(^max) 



W 2 (k n 



(51) 



The assumptions of Eq. (49), (50) and (51) are illustrated in 
Fig. 1, in which we show the simulated power spectrum for 
the Planck High Frequency Instrument (HFI) 100 GHz fre- 
quency channel. Eq. (50) and Eq. (51) lead to the following 
relationship between fc max and k m in parameter, 



-S(&max)^(^><min) 



c pi: 



C- k y(k min )' 



which leads to 



9 2 
h ~ — In 

^max — q2 



mm / 



(52) 



(53) 




Figure 3. The power spectrum of the foreground-cleaned map by 
Tegmark, de Olivcira-Costa & Hamilton (2003) and the WMAP 
Q-band map. The shaded areas are the filtering ranges of the 
ATHF for point source detection: (£ m i n ,&nax) = (200,400) for 
the FCM and (400, 900) for the WMAP Q-band map. 



nATHF 2 \ ^ a 2 

R - g y J Si- 



»ic^(fc min ) • 



(54) 



We can further solve Eq. (47) by inserting Eq. (27), which 
then becomes (with the integral part of C pix < C sky W 2 ) 



f 

Jk„ 



dk k C sky {k)W\k) 



C ^(/Cmin) 



2<9 2 



(55) 



If the C sky (k) = Ak n , Eq. (55) can be approximated and 
we reach 



kit 



II I 1 u i I- 



(56) 



The ratio 7? can then be obtained from Eq.(44) 



Note that we use the simple trapezoidal rule for the integra- 
tion, which works better when the index n is smaller. 



4 APPLICATION OF THE ATHF ON 
SIMULATIONS AND WMAP DATA 

4.1 Simulations of the Planck HFI 100 GHz 
channel 

We carry out the simulations of the Planck HFI 100 GHz 
channel for point source extraction by the TOF and the 
ATHF. The simulation area is 25.6 deg 2 , and the pixel 
size is 3.0 arcmin. The ctcmb and the <r p i x are 4.03 x 10~ 5 
and 6.07 x 10 -6 , respectively. Point sources with amplitudes 
above 2ctcmb at this channel would be flushed out by both 
(symmetric) filters with the 5 a criterion (see the 5th col- 
umn of Table 3 in Chiang et al. (2002b)), so we add ran- 
domly 50 point sources whose amplitudes are 1.5ctcmb to 
check the efficiency. Obviously the number of point sources 
is too large for this frequency channel, but it serves as the 
demonstration of the efficiency of extraction. We specifically 
use an elliptical beam shape in the simulation with elliptic- 
ity ratio 0+/6- = 1.3 (Burigana et al. 1998). We desig- 
nate the elliptical beam size such that 9 2 _ + 9% = 29q where 
0o = FWHM/V8 In 2. For the HFI 100 GHz channel the 
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Figure 4. The adaptive top-hat filtered maps from the foreground-cleaned map (top) by Tegmark, de Oliveira-Costa & Hamilton (2003) 
and the WMAP Q-band map (bottom). The filtering ranges for both maps are shown in Fig. 3. In the top panel, the circled are peaks 
with amplitudes above 5<Tf after filtering, which are the point source residues from the component separation. One interesting feature 
is that not all the residues we recover are extragalactic point sources. The dash circles are not peaks, but dips. In the bottom panel we 
retain the positions of the circles for comparison, from which those peaks are point sources. This is to prove that most of the peaks we 
retrieve from the FCM are point source residues. 
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FWHM is 10.7 arcmin. In order to simulate more realis- 
tic situations, we put this beam in rotation while scanning 
across the simulation area in order to investigate how the 
elliptical beam could affect the efficiency of point source ex- 
traction for the TOF and the ATHF. 

The period of rotation of the beam is 2n along both 
sides of the simulation square. Although the rotation period 
is too large for such size of map, it nevertheless can elucidate 
the effect of elliptical beam shape regarding point source 
extraction. 

Without the exact information of the in-flight beam size 
and orientation, the optimal form of the TOF is by insert- 
ing a circular beam function into Eq. (41), where the win- 
dow function Wl = H\ = exp(-fc 2 6» 2 ). In Fig. 2 we plot 
the point source extraction rate from the TOF for the 50 
point sources. The plot is the extraction percentage of point 
sources against the 6 inserted to the TOF (in terms of 100 
GHz channel 9q). We would like to examine, for in-flight el- 
liptical beam shapes, how the supplied beam function for the 
TOF can affect the extraction rate. In the dotted curve, the 
power spectrum is precisely known. The short-dash curve is 
the case when the slope of the CMB power spectrum has 5 
per cent more tilt, whereas the long-dash line 5 per cent less 
tilt. The horizontal line marks the ATHF extraction per- 
centage. With the exact power spectrum, the TOF is able 
to detect 85 per cent of the point sources with amplitude 
1.5o"cmb when the supplied 9 is chosen in between 9+ and 
6- . For the case of 5 per cent more tilt in power spectrum, 
the supplied beam function will need adjustment to a smaller 
9 value to ensure the TOF to have a maximal resonance, 
and adjustment to a larger 9 for the case of 5 per cent less 
tilt. On the other hand, the ATHF can reach 74 per cent 
of extraction rate without any information. From Fig. 2 it 
is not surprising that the 'safe bet' of the beam function is 
9/9o — 1 for the uncertainties in the power spectrum. How- 
ever, as the cleaning of foreground contamination precedes 
the determination of the angular power spectrum, and dur- 
ing the whole sky scan of the Planck mission, the possible 
degradation effect of the mirrors could change the beam size 
and shape, the efficiency of the TOF will be hampered. 



4.2 The ATHF on the WMAP derived map 

In this subsection we apply the ATHF to the derived map 
from the WMAP 1-year data~t\ The derived maps are pro- 
duced by Tegmark, de Oliveira-Costa & Hamilton (2003) 
performing an independent component separation analysis 
from the WMAP team. They are the foreground-cleaned 
map (FCM) and the Wiener-filtered map (WFM). The FCM 
was tested to be non-Gaussian (Chiang, Naselsky, Verkho- 
danov & Way 2003), partly due to galactic emission. Here 
we investigate the possible point source residues from their 
component separation. 

We use the ATHF on the FCM as it is the map with 
scientific significance. In Fig. 3 we show the power spectra 
of the FCM and the WMAP Q-band map with the corre- 
sponding filtering ranges (the shaded area) of the ATHF: 
(l mi „,f ma x) = (200,400) for the FCM and (400,900) for the 

i http: / /lambda.gsfc. nasa.gov/product/wmap/ 



WMAP Q-band map. It is easy to see that the CMB (plus 
foreground) power spectra (convolved with the beam) and 
the noise level. Note that the presentation of the angular 
power spectrum in Fig. 3 is the standard format, i.e. Ct 
with the factor 1(1 + l)/47r, which is different from Fig. 1 
used for analysis in the previous Section. We therefore apply 
the ATHF with filtering ranges covering the conjunction of 
the CMB and noise power curves. The choice of the filtering 
range follows the rule of the thumb described in Chiang et 
al. (2002b): the £ m i n to exclude most of the CMB power, 
the €max the noise power, and both to keep the part being 
convolved by the beam. 

Fig. 4 shows the filtered maps of the FCM (top panel) 
and the WMAP Q-band map (bottom). The circled in the 
top panel are filtered peaks with amplitudes above 5<7f . One 
interesting feature in Fig. 4 is that not all the residues we 
recover are extragalactic point sources. The dash circles are 
not peaks, but dips. In the bottom panel we retain the posi- 
tions of the circles for comparison, from which those peaks 
are point sources. This is to prove that most of the peaks we 
retrieve from the FCM are point source residues. We do not 
intend to recover all residues from one single iteration of the 
ATHF but simply demonstrate the usefulness and quickness 
of the ATHF on the real data. Application of the ATHF on 
the WMAP 5 different channel maps will be on a separate 
paper. 

One important issue related to elliptical beam- 
convolved point source extraction is the removal of point 
sources. Normally, circular Gaussian profile is used for clean- 
ing the elliptical beam- convolved point sources from the map 
after detection of the point source position. Due to the el- 
lipticity of the beam and pixel noise, there are considerable 
residues after the cleaning by circular Gaussian profile. To 
eliminate the residual peaks, we apply a simple algorithm 
that is modified from the so-called hybrid median filter used 
in signal processing. The hybrid median filter is useful in 
noise reduction and in particular in elimination of shot noise. 
This non-linear filter replaces the targeting pixel with aver- 
age of the neighboring pixels from the same row, column 
and two diagonals. In the modified cascading version, the 
targeted region (3x3 pixels) centered at the residual peak 
is replaced with the following algorithm (see Fig. 5). The 
central pixel is replaced by the median of the 8 pixels from 
the same row, column, and two diagonals right outside the 
targeted region. The four pixels at the four corners in the 
region can then be decided by the median of their 4 cor- 
ner pixels, which include the already-replaced central pixel. 
The rest 4 pixels can be replaced by the median of their 
neighbouring pixels from the same row and column. 



5 CONCLUSIONS 

We have discussed analytically the three main linear filters 
for point source extraction of the Planck mission. The TOF 
is optimal in terms of the gain factor. The so-called optimal 
pseudo-filter is at its best asymptotic at the far ends to 
the TOF under certain circumstance. Both of these filters 
require the experiment parameters such as the beam shape 
and size. 

Note that the calculation is based on the simple geomet- 
rical model of the beam. As we mention in the beginning, 
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Figure 5. The cascading hybrid median filtering scheme. For 
a residual peak located at from the cleaning by circular 

Gaussian profile, we can apply the cascading hybrid median filter 
to replace the inner 3x3 pixels. First of all, we replace the central 
pixel by the median of the 8 pixels from the same row, 

column, and the diagonals outside the targeted region, i.e., by 
the median of pixels (i-2,j), (i-2,j + 2), (i,j + 2), (i + 2,j + 2), 
(i + 2,j), (i + 2,j - 2), - 2) and (i - 2, j - 2). Once the 
central pixel is decided, we can in turn replace the four corners 
by the median of its diagonals, e.g., the pixel (i + l,j + 1) can be 
replaced by the median of + 2), (i + 2,j + 2), (i + 2,j) and 
the replaced (i, j). Then the last 4 pixels can be replaced by their 
neighbouring pixels of the same row and column, e.g., (i + 1, j) 
by the median of (i + 1, j + 1), (i + 2,j), (i + l,j — 1) and the 
replaced (i, j). 



the point source is the reflection of, among others, the beam 
shape properties. The possible antenna beam degradation 
would be a problem for the calibration of the in-flight an- 
tenna beam shape (Naselsky et al. 2002) , causing the change 
of the inclination and the ellipticity ratio of the beam during 
the mission. 

There is one comment on the extraction of point sources 
from the time-ordered scans when the inclination of the 
beam, the issue of pointing and the noise properties are most 
crucial: the inclination angle and the pointing will decide the 
point source FWHM on the time-ordered data, for which the 
ATHF is particularly useful. 

Note also that we can easily generalize the multi- 
frequency method for point source extraction (Naselsky, 
Noviko & Silk 2002) by applying the ATHF for all Planck 
frequency channels. 
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